SoSe2021

Folienübersicht

Kontinuierliche und kategoriale X: Kovarianzanalyse

Kovarianzanalyse oder ANalysis of COVAriance (ANCOVA)

Ein Hybrid zwischen Regression und ANOVA

  • Einige betrachten die ANCOVA als ein “ANOVA”-Modell mit einer kontinuierlichen X Variable (Kovariate) → Fokus liegt auf den Faktorstufen, welche für die Kovariate angepasst wurden (z. B. Quinn & Keough, 2002).
  • Andere betrachten es als ein “Regressions”-Modell mit einem kategorialen Prädiktor → Fokus liegt auf der Kovariate, welche für den Faktor angepasst wurde (z. B. Crawley 2005,2007).
  • Das typische “Maximal”-Modell beinhaltet die Schätzung eines Achsenabschnitts und einer Steigung (Regressionsteil) für jede Stufe der kategorialen Variable(n) (ANOVA-Teil).

Hybrid zwischen Regression und ANOVA

Gewichtsbeispiel 1

Aus der Regressionsperspektive

\[Weight_{ij}=\alpha_i+\beta_i*Age_j + \epsilon_{ij}\]

Gleichung nach Geschlecht getrennt

\[Weight_{female,j}=\alpha_{female}+\beta_{female}*Age_j + \epsilon_{female,j}\] \[Weight_{male,j}=\alpha_{male}+\beta_{male}*Age_j + \epsilon_{male,j}\]

Hybrid zwischen Regression und ANOVA

Gewichtsbeispiel 2

Mögliche Modelltypen

Gewichtsbeispiel 1

Mögliche Modelltypen

Gewichtsbeispiel 2

Mögliche Modelltypen

Gewichtsbeispiel 3

Signifikanztests

Achsenabschnitte \(\alpha_i\)

  • Nullhypothese \(H_0: \alpha_1 = \alpha_2 = ... \alpha_i\)
  • Bei 2 Achsenabschnitten: t-Test
  • Bei >2 Achsenabschnitten: ANOVA (F-Test)

Steigungen \(\beta_i\)

  • Nullhypothese \(H_0: \beta_1 = \beta_2 = ... \beta_i\)
  • Bei 2 Steigungen: t-Test
  • Bei >2 Steigungen: ANOVA (F-Test)

Beispiel aus DS1: Königslachse

ChinookArg

Logarithmierte Gewichtsunterschiede und Längen-Gewichtsbeziehung

Längen und Gewichte für Königslachse von drei Orten in Argentinien (Daten sind im FSA Paket enthalten):

data("ChinookArg", package = "FSA")
str(ChinookArg)
'data.frame':   112 obs. of  3 variables:
 $ tl : num  120 115 111 110 110 ...
 $ w  : num  17.9 17.2 16.8 15.8 14.3 13.8 12.8 11.7 12.8 14.8 ...
 $ loc: Factor w/ 3 levels "Argentina","Petrohue",..: 1 1 1 1 1 1 1 1 1 1 ...

Da für die meisten Arten die Beziehung zwischen dem Gewicht eines Tieres und seiner Länge einer 2-Parameter Powerfunktion (\(W = aL^{b}\)) folgt, wollen wir die Variablen log-transformieren:

ChinookArg$tl_log <- log(ChinookArg$tl)
ChinookArg$w_log <- log(ChinookArg$w)

Königslachse

Datenexploration

Wie ist die Längen-Gewichtsbeziehung und gibt es Ortsunterschiede?

Königslachse

Maximalmodell

full_mod <- lm(w_log ~ loc + tl_log + loc:tl_log, ChinookArg) # oder: ~ loc * tl_log
summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.58273 -0.18471 -0.00186  0.13088  1.63620 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -6.6750     1.5904  -4.197 5.64e-05 ***
locPetrohue         -2.3957     3.1494  -0.761    0.449    
locPuyehue          -2.0696     1.6868  -1.227    0.223    
tl_log               1.9836     0.3530   5.619 1.56e-07 ***
locPetrohue:tl_log   0.4795     0.6928   0.692    0.490    
locPuyehue:tl_log    0.3624     0.3793   0.955    0.342    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8924 
F-statistic:   185 on 5 and 106 DF,  p-value: < 2.2e-16

Kleiner Exkurs:
Dummy-Variablen in R

Wenn erklärende Variablen kategorial sind

Falls x kategorial ist, macht es nicht viel Sinn das Model \(y = a + bx\) zu benutzen, da \(x\) nicht mit \(b\) multipliziert werden kann. Aber R hat eine Hilfskonstruktion:

  • Die einzelnen Stufen von kategoriale Variablen werden in sog. Dummy-Variablen jeweils konvertiert.
  • Jede Dummy-Variable ist als 0 (FALSE) oder 1 (TRUE) kodiert.
  • Anzahl der Dummy-Variablen = Anzahl der Gruppen minus 1
  • Alle linearen Modelle passen kategoriale erklärende Variablen mit Hilfe von Dummy-Variablen an.

Demonstration mit model.matrix()

Die Daten

df
      y length
1   6.7      M
2  16.6      M
3  12.8      M
4  16.1      M
5   4.2      L
6   1.3      S
7  18.6      L
8   5.0      S
9   8.4      S
10  2.0      S

Die Dummy-Variablen

modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
   `(Intercept)` lengthM lengthS
           <dbl>   <dbl>   <dbl>
 1             1       1       0
 2             1       1       0
 3             1       1       0
 4             1       1       0
 5             1       0       0
 6             1       0       1
 7             1       0       0
 8             1       0       1
 9             1       0       1
10             1       0       1

Wo ist die Längen-Klasse L geblieben?

Demonstration mit model.matrix()

Die Daten

df
      y length
1   6.7      M
2  16.6      M
3  12.8      M
4  16.1      M
5   4.2      L
6   1.3      S
7  18.6      L
8   5.0      S
9   8.4      S
10  2.0      S

Die Dummy-Variablen

modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
   `(Intercept)` lengthM lengthS
           <dbl>   <dbl>   <dbl>
 1             1       1       0
 2             1       1       0
 3             1       1       0
 4             1       1       0
 5             1       0       0
 6             1       0       1
 7             1       0       0
 8             1       0       1
 9             1       0       1
10             1       0       1

L ist durch den Achsenschnittpunkt (erste Spalte) dargestellt!

Dummy-Variablen am Beispiel der Königslachse

model.matrix(w_log ~ loc, ChinookArg) %>% 
    as.data.frame(.) %>%
  # Hinzufügen der Original loc Variable
  mutate(sampl_loc = ChinookArg$loc) %>%
  head()
  (Intercept) locPetrohue locPuyehue sampl_loc
1           1           0          0 Argentina
2           1           0          0 Argentina
3           1           0          0 Argentina
4           1           0          0 Argentina
5           1           0          0 Argentina
6           1           0          0 Argentina

Wir können sehen, dass die ersten 5 Lachse in ‘Argentinien’ gefangen wurden und dass die Beobachtungen mit einer 1 für die ‘intercept’ Variable (erste Faktorstufe → hier Argentinien) kodiert und 0 für die Variable locPetrohue und locPuyehue wurden.

Interpretation der Koeffizienten

1

Wie würdest Du diese Koeffizienten interpretieren (wenn also nur der Fangort die erklärende Variable ist)?

chin_mod <- lm(w_log ~ loc, ChinookArg)
coef(chin_mod)
(Intercept) locPetrohue  locPuyehue 
 2.25605528 -0.09844412 -1.52993775 

Wie lautet die Gleichung für den Fangort?

Interpretation der Koeffizienten

2

Setze die berechneten Koeffizienten in die Regressionsgleichung mit Dummy-Variablen ein:

Zurück zum Maximalmodell mit Interaktion

1

summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.58273 -0.18471 -0.00186  0.13088  1.63620 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -6.6750     1.5904  -4.197 5.64e-05 ***
locPetrohue         -2.3957     3.1494  -0.761    0.449    
locPuyehue          -2.0696     1.6868  -1.227    0.223    
tl_log               1.9836     0.3530   5.619 1.56e-07 ***
locPetrohue:tl_log   0.4795     0.6928   0.692    0.490    
locPuyehue:tl_log    0.3624     0.3793   0.955    0.342    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8924 
F-statistic:   185 on 5 and 106 DF,  p-value: < 2.2e-16

Zurück zum Maximalmodell mit Interaktion

2

Grundgleichung

\[w.log_{ij} = \alpha_{i} + \beta_{i}*tl.log_{j} + \epsilon_{ij}\]

wobei i = Fangort, j = jeder einzelne Lachs


Ortspezifische Gleichungen

\[w.log_{A,j} = \alpha_{A} + \beta_{A}*tl.log_{j} + \epsilon_{A,j}\] \[w.log_{Pe,j} = \alpha_{Pe} + \beta_{Pe}*tl.log_{j} + \epsilon_{Pe,j}\] \[w.log_{Pu,j} = \alpha_{Pu} + \beta_{Pu}*tl.log_{j} + \epsilon_{Pu,j}\]

Berechnung der mittleren Gewichte pro Fangort

Gleichung mit Dummy Variablen

\[w.log = a_A + a_{Pe}*locPetrohue + a_{Pu}*locPuyehue +\\ b_A*tl.log+ b_{Pe}*locPetrohue*tl.log + b_{Pu}*locPuyehue*tl.log\]

Einsetzen der Koeffizienten aus R in die Gleichungen

coef(full_mod)
       (Intercept)        locPetrohue         locPuyehue             tl_log 
        -6.6749719         -2.3957331         -2.0695530          1.9836016 
locPetrohue:tl_log  locPuyehue:tl_log 
         0.4794543          0.3624014 

Argentina: \(w.log_{A} = -6.67+1.98*tl_log\)

Petrohuye: \(w.log_{Pe}= (-6.67+-2.40)+(1.98+0.48*1)*tl.log = -9.07+2.46*tl.log\)

Puyehuye: \(w.log_{Pu} = (-6.67+-2.07)+(1.98+0.36*1)*tl.log = -8.74+2.34*tl.log\)

Systematische Modellauswahl

1

Anstelle der ‘all subsets selection’ (siehe Bsp. airquality), wollen wir hier etwas systematischer vorgehen und die sog. ‘backward selection’ vornehmen.

‘backward selection’

  1. Starte mit dem Maximalmodell, welches alle X Variablen und mögliche Interaktionen enthält, und überprüfe die Modellannahmen anhand dieses vollen Modells.
  2. Erstelle nun mehrere reduzierte Modelle, bei der jeweils 1 Interaktionsterm oder ein Hauptterm (der nicht in einem Interaktionsterm vorkommt) entfernt wurde.
  3. Vergleiche das volle mit den reduzierten Modellen mittels z.B. AIC und wähle das Modell, welches am meisten erklären kann bzw. den niedrigsten AIC hat: dieses stellt nun das neue ‘volle’ Modell dar.
  4. Erstelle erneut reduzierte Modelle und vergleiche mittels AIC.
  5. Wiederhole diesen Prozess bis das ‘volle’ Modell den niedrigsten AIC hat.

Systematische Modellauswahl

2

Alternativen

  • ‘forward selection’: hier wird genau umgekehrt verfahren, es wird mit dem Nullmodell gestartet und nacheinander eine Variable hinzugefügt.
  • ‘stepwise forward/backward selection’:
    • Kombiniert beide Ansätze.

R Funktionen

  • step() → führt eine vollständige Selektion durch, bis das Modell mit dem niedrigsten AIC gefunden wird.
    • Mit dem Argument method kann eingestellt werden, ob die Richtung ‘forward’, ‘backward’ oder beides (‘both’) sein soll.
    • Hinweis: Die Funktion wendet nicht die Differenzregel (von mind. 2 AIC Werten) für komplexere Modelle an, d.h. es muss hier evt. manuell weiter selektiert werden!
  • drop1 → erstellt reduzierte Modelle des übergebenen vollen Modells und gibt eine Vergleichstabelle mit den AICs zurück.

Königslachse

Prüfen der Annahmen am vollen Modell 1

par(mfrow = c(2,2))
plot(full_mod)

Königslachse

Prüfen der Annahmen am vollen Modell 2

Zusätzliche Plots: Residuen vs. jede X Variable

ChinookArg$residuals <- resid(full_mod)
par(mfrow = c(1,2), mar = c(4,4,1,1))
boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2)
plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)

Königslachse

Modellreduktion - Manuell

Wenn Interaktionsterme im Modell vorkommen, müssen auch die enthaltenen Variablen als Hauptterme vorkommen. Es darf also nicht tl_log oder loc im reduzierten Modell fehlen. Nur die Interaktion selbst kann aktuell rausgenommen werden:

Runde 1

red_mod <- update(full_mod, .~.-loc:tl_log)
AIC(full_mod, red_mod)
         df      AIC
full_mod  7 70.53732
red_mod   5 67.57480
  • Das reduzierte Modell hat den niedrigeren AIC, es kann also der Interaktionsterm entfernt werden.

Runde 2

# Das reduzierte Modell wird das volle:
full_mod <- red_mod 

# Einzelne X entfernen
red_mod1 <- update(full_mod, .~.-loc)
red_mod2 <- update(full_mod, .~.-tl_log)

AIC(full_mod, red_mod1, red_mod2)
         df       AIC
full_mod  5  67.57480
red_mod1  3  87.69126
red_mod2  4 224.06339
  • Das volle Modell hat den besten AIC, daher wird nichts weiter entfernt.

Königslachse

Modellreduktion - automatisch

Built-in Funktion step()

step(full_mod, method = "backward")
Start:  AIC=-249.3
w_log ~ loc + tl_log + loc:tl_log

             Df Sum of Sq    RSS     AIC
- loc:tl_log  2    0.1011 10.965 -252.27
<none>                    10.864 -249.31

Step:  AIC=-252.27
w_log ~ loc + tl_log

         Df Sum of Sq    RSS      AIC
<none>                10.965 -252.267
- loc     2     2.634 13.599 -232.151
- tl_log  1    34.175 45.140  -95.779
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)

Coefficients:
(Intercept)  locPetrohue   locPuyehue       tl_log  
    -8.1217      -0.2281      -0.4570       2.3049  

  • Zum Schluss wird das ‘optimale’ Modell zusammengefasst.
  • Auch hier wird das Modell mit beiden Variablen ohne Interaktion als ‘bestes’ Modell gewählt.

Königslachse

Validierung des finalen Modells 1

best_mod <- lm(w_log ~ loc + tl_log, ChinookArg)
par(mfrow = c(2,2))
plot(best_mod)

Königslachse

Validierung des finalen Modells 2

Zusätzliche Plots: Residuen vs. jede X Variable

ChinookArg$residuals <- resid(best_mod)
par(mfrow = c(1,2), mar = c(4,4,1,1))
boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2)
plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)

Königslachse

Ergebnis 1 - Regressionkomponente

summary(best_mod)
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.60633 -0.19351 -0.00076  0.14352  1.61286 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.12169    0.56827 -14.292  < 2e-16 ***
locPetrohue -0.22813    0.08013  -2.847  0.00528 ** 
locPuyehue  -0.45699    0.09231  -4.951 2.74e-06 ***
tl_log       2.30492    0.12563  18.347  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3186 on 108 degrees of freedom
Multiple R-squared:  0.8962,    Adjusted R-squared:  0.8934 
F-statistic:   311 on 3 and 108 DF,  p-value: < 2.2e-16

Königslachse

p-Werte für kategoriale Variablen

Problem der t-Tests und p-Werte

# Extraktion der Koeffizientenschätzung mit tidy() aus dem broom Paket
broom::tidy(best_mod)
# A tibble: 4 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   -8.12     0.568     -14.3  1.19e-26
2 locPetrohue   -0.228    0.0801     -2.85 5.28e- 3
3 locPuyehue    -0.457    0.0923     -4.95 2.74e- 6
4 tl_log         2.30     0.126      18.3  5.72e-35
  • Diese basieren auf dem Vergleich der einzelnen Faktorstufen mit der ‘Intercept’, also der Baseline-Faktorstufe.
  • Möglicherweise unterscheiden sich diese nicht zur Baseline, aber untereinander!

Königslachse

Ergebnis 2 - ANOVA-Komponente

Um die ANOVA Tabelle zu erhalten, wo alle Faktorstufen miteinander verglichen werden, nutze die Funktion anova() für dein lm Objekt:

anova(best_mod)
Analysis of Variance Table

Response: w_log
           Df Sum Sq Mean Sq F value    Pr(>F)    
loc         2 60.542  30.271  298.16 < 2.2e-16 ***
tl_log      1 34.175  34.175  336.61 < 2.2e-16 ***
Residuals 108 10.965   0.102                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

→ Dieser p-Wert für loc sollte im Bericht stehen, nicht die aus der summary Funktion!

Königlachse

Gleichung

coef(best_mod)
(Intercept) locPetrohue  locPuyehue      tl_log 
 -8.1216871  -0.2281264  -0.4569901   2.3049203 

Argentina: \(w.log_{A} = -8.12+2.30*tl.log\)

Petrohuye: \(w.log_{Pe}= (-8.12+-0.23)+2.30*tl.log = -8.35+2.30*tl.log\)

Puyehuye: \(w.log_{Pu} = (-8.12+-0.46)+2.30*tl.log = -8.58+2.30*tl.log\)

Königslachse

Visualisierung

ChinookArg %>% data_grid(
  tl_log = seq_range(tl_log, 20),
  loc = factor(levels(loc))
) %>%
  add_predictions(best_mod) %>%
  ggplot(aes(x = tl_log, group = loc, 
    colour = loc)) + 
  geom_line(aes(y = pred)) +
  geom_point(data = ChinookArg, 
    mapping = aes(y = w_log)) +
  ylab("w_log") +
  scale_colour_manual(values = 
      c("forestgreen", "orange1", 
        "deeppink4")) +
  theme_classic()

Ablaufprotokoll

Ablaufprotokoll

Beispiel für ‘backward selection’

  1. Datenexploration:
    • Ausreißer in einer der X Variablen? Variable entfernen oder transformieren.
    • Ausreißer in Y? Behalten und hoffen, dass diese nicht zu einflussreich sind.
    • Y nicht normalverteilt? Beibehalten und hoffen, dass die Residuen normalverteilt sind.
    • Zusammenhang nicht linear? Variablen transformieren (entweder X, oder Y oder beide).
    • X kollinear? Variable mit höchstem VIF-Wert entfernen oder nach biologischem Wissen entscheiden.
    • Interaktionen? Evtl. explizit mit modellieren.
  2. Starte mit dem vollen Modell und beziehe Interaktionen nur mit ein, wenn Du vorher weißt, dass sie eine Rolle spielen könnten → nicht auf numerischen Output schauen → erst Annahmen überprüfen: wenn starke Heterogenität in den Residuen ggf. (weitere) Transformation.
  3. Führe die Modellauswahl durch, d.h. entferne schrittweise einzelne Variablen bis der AIC minimal ist.
  4. Annahmen des “optimalen” Modells prüfen.
  5. Wenn alles ok, weiter zu 8.
  6. Wenn nicht:
    • Residuen nicht normal? Keine Sorge, Homogenität ist wichtiger.
    • Irgendwelche einflussreichen Beobachtungen? Variable(n) entfernen, Beobachtungen (ganze Zeilen) entfernen oder transformieren.
    • Heterogenität, d.h. Muster im Streudiagramm Residuen vs. Fit?
    • Abhängigkeit, d. h. Muster im Streudiagramm Residuen vs. jedes X (auch die im optimalen Modell ausgeschlossenen)?
  7. Passe das Modell erneut an und überprüfe wieder die Annahmen.
  8. Wenn diese nun erfüllt sind, betrachte das numerische Ergebnis: wie viel Variabilität wird erklärt? Ist der Einfluss von X signifikant? Wie stark ist der Einfluss?
  9. Stelle das Ergebnis in einer Tabelle und in einem Streudiagramm dar (einschließlich der Regressionslinie).

Muster in den Residuen und mögliche Lösungen

  1. Keine Homogenität & kein klares Muster in den Residuen gegen jedes X:
    • Transformiere die Antwortvariable.
    • Füge Interaktionen hinzu.
    • Füge weitere erklärende Variablen hinzu.
    • Anwendung generalisierter linearer Modelle (GLM) mit z.B. Poisson-Verteilung wenn die Daten Zählungen sind.
  2. Keine Homogenität & klares Muster in den Residuen gegen jedes X:
    • Füge Interaktionen hinzu.
    • Wenn die erklärende Variable, bei der ein Muster zu erkennen ist, aus dem endgültigen Modell ausgeschlossen wurde, wieder hinzufügen.
    • Füge nicht-lineare Terme der erklärenden Variable hinzu (z. B. quadratische Terme).
    • Anwendung generalisierter additiver Modelle (GAM).
  3. Homogenität gegeben, aber klares Muster in den Residuen gegen jedes X:
    • Transformiere die entsprechenden erklärenden Variablen
    • Wenn die erklärende Variable, bei der ein Muster zu erkennen ist, aus dem endgültigen Modell ausgeschlossen wurde, wieder hinzufügen.
    • GAM anwenden.

Fragen?